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Introduction 

Sensitivity analysis is emerging as a fruitful area of engineering research. The rea- 
son for this interest is the recognition of the variety of uses for sensitivity derivatives. In 
its early stages, sensitivity analysis found its predominant use in assessing the effect of 
varying parameters in mathematical models of control systems; see, for example, the texts 
of Toinovic, 1 Brayton and Spence, 2 Prank, 3 and Radanovic 4 for discussions of the early 
development of sensitivity theory. Interest in optimal control in the early 1960s (see, for 
example, Ref. 5) and automated structural optimization (see for example, Ref. 6) led to the 
use of gradient-based mathematical programming methods in which derivatives were used 
to find search directions toward optimum solutions. More recently, there has been strong 
interest in promoting systematic structural optimization as a useful tool for the practicing 
structural design engineer on large problems — a process still under way. Early attempts 
to use formal optimization for large structural systems resulted in excessively long and 
expensive computer runs. Examination of the optimization procedures indicated that the 
predominant contributor to the cost and time was the calculation of derivatives. As a con- 
sequence, emerging interest in sensitivity analysis has emphasized efficient computational 
procedures. In addition, researchers have developed and applied sensitivity analysis for 
approximate analysis, analytical model improvement, and assessment of design trends — so 
that structural sensitivity analysis has become more than a utility for optimization and is 
a versatile design tool in its own right. Most recently, researchers in disciplines such as 
physiology, 7 thermodynamics, 8 physical chemistry, 9 and aerodynamics, 10-12 have been us- 
ing sensitivity methodology to assess the effects of parameter variations in their analytical 
models and to create designs insensitive to parameter variation. 13,14 

Derivatives of structural response can be calculated analytically at three stages. We 
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can differentiatie the continuum equations defining the response of the structure. We can 
differentiate the equations obtained when the continuum equations are discretized which 
is the topic of the present chapter. Finally, we can differentiate directly the computer 
program used to solve the structural response, such as a finite element computer program. 
This third approach is not discussed in this textbook, but the interested reader is referred 
to Refs 15 and 16. Analytical derivative calculations typically entail a substantial effort of 
analysis and software development. In many cases it is better to use derivatives obtained 
from a finite difference approximation. This chapter therefore starts with the discussion 
of the calculation of derivatives by finite differences. 

Finite Difference Sensitivities 
Truncation and Condition Errors 

The simplest finite difference approximation is the first-order forward-difference ap- 
proximation. Given a function g(v ) of a design variable v, the forward-difference approxi- 
mation Ag/Av to the derivative dg/dv is given as 


A 9 = g(t> + An) - ff(u) m 

Av Av K ’ 

Another commonly used finite-difference approximation is the second-order central-difference 
approximation 


A g _ g(u + A v) - g(v - Av ) 

Av 2Av 

It is also possible to employ higher-order finite-difference approximations, but they are 
rarely used in structural optimization applications because of the associated high compu- 
tational cost. If we need to find the derivatives of the structural response with respect to n 
design variables the forward-difference approximation requires n additional analyses, the 
central-difference approximation 2 n additional analyses, and higher order approximations 
are even more expensive. 

The key to the selection of the approximation and the step size Ax is an estimate of 
the required accuracy. This topic is discussed in Ref. 17, and is summarized next. 

Whenever finite-difference formulae are used to approximate derivatives, there are 
two sources of error: truncation and condition errors. The truncation error e T (Av) is a 
result of the neglected terms in the Taylor series expansion of the perturbed function. For 
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example, the Taylor series expansion of g(v + Av) can be written as 


g{v + Av ) - g(v ) + Av^(v) + — (v 4- CAn), 0<C<1 (3) 

From Eq. (3) it follows that the truncation error for the forward-difference approximation 
is 

e T (Av) = + (Av) 0<C<1 (4) 

Similarly, by including one more term in the Taylor series expansion we get that the 
truncation error for the central difference approximation is 


e r (Av) = — (u + CAv) 0<C<1 (5) 

The condition error is the difference between the numerical evaluation of the function and 
its exact value. One contribution to the condition error is round-off error in calculating 
dg/dv from the original and perturbed values of g. This contribution is comparatively 
small for most computers unless At> is extremely small. However if g(x) is computed by 
a lengthy or ill-conditioned numerical process, the round-off contribution to the condition 
error can be substantial. Additionally, condition errors may result if g(x) is calculated 
by an iterative process which is terminated early. If we have a bound e g on the absolute 
error in the computed function g , we can estimate the condition error. For example, for 
the forward-difference approximation the condition error ec(Av) is (very!) conservatively 
estimated from Eq. (1) as 


ec( Av) = -^e 9 (G) 

Equations (4) and (6) present us with the so called “step-size dilemma.” If we select the 
step size to be small, so as to reduce the truncation error, we may have an excessive 
condition error. In some cases there may not be any step size which results in acceptable 
error! 

A bound e on the total error, the sum of the truncation and condition errors, for the 
forward-difference approximation is obtained from Eqs. (4) and (6) as 


e = 


At), , 2 

T w + s; e « 


(7) 
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where Sb is a bound on the second derivative in the interval [v,v + Av]. When e g and Sb 
are available it is possible to calculate an optimum step-size that minimizes e as 



( 8 ) 


Procedures for estimating s& and e u are given in Ref 17 and 18. 

Iteratively solved Problems 

Condition errors can become important when iterative methods are used for perform- 
ing some of the calculations. Consider a simple example of a single displacement component 
u which is obtained by solving a nonlinear algebraic equation which depends on one design 
variable v 

f(v y u) = 0 (9) 

The solution of Eq. (9) is obtained by an iterative process which starts with some initial 
guess of u and terminates when the iterant u is estimated to be within some tolerance e of 
the exact u (Note that e is a bound on the condition error in tt). To calculate the derivative 
du/dv assume that we use the forward-difference approximation. That is, we perturb v by 
Av and solve Eq. (9) for 

f(v + An, ua) = 0 (10) 

The iterative solution of Eq.(10) yields an approximation and then du/dv is approxi- 
mated as 

ftli 'Tl A — 7/ 

( 11 ) 


du 

dv 


u A 


Av 


To start the iterative process for obtaining ua, two initial guesses come to mind. The first 
is to start with the same initial guess that was used to solve for u. If the convergence of 
the iterative process is monotonic there is a good chance that when we use Eq. (11) the 
errors in u and ua will almost cancel out, and we will get a very small condition error. 
The other logical initial guess for it a is u. This initial guess is known to be good because 
Ax is typically small, and so we may get fast convergence. Unfortunately, this time we 
cannot expect the condition errors to cancel. As we iterate on the original error (the 
difference between u and u) will be reduced at the same time that the change due to Ax 
is taking effect (consider, for example, what happens if Ax is set to zero, or an extremely 
small number). 

Reference 19 suggests a strategy which allows us to start the iteration for ua from u 
without worry of excessive condition errors. The approach is to pretend that u is the exact 
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rather than approximate solution by changing the problem that we want to solve. Indeed, 
u is the exact solution of 

f(v, u) - f(v, u) = 0 (12) 

which is only slightly different from our original problem (because f(v, u) is almost zero). 
We now find the derivative du/dv from Eq.(12), by obtaining u& as the solution of 

f(v + Aw, u A ) - f(v, u) = 0 (13) 

Because u is the exact solution of this equation for An = 0 the iterative process will only 
reflect the effect of Aw, and we will obtain a good approximation from Eq. (11). 

Because of the high cost and the accuracy problems associated with finite-difference 
derivatives, there has been much effort into developing analytical derivative approxima- 
tions. The rest of this chapter is devoted to such analytical expressions for sensitivity 
derivatives. 


Sensitivity of Static Response 
First Derivatives of Linear Response 

This section of the paper focuses on the calculation of first derivatives of static linear 
structural response (displacements and stresses) computed from finite element models. 
The governing equation for displacement is 

. KU = F (14) 

where K is the symmetric stiffness matrix of order n x n, U the vector of displacements, 
and F the vector of applied forces. Both K and F are, in general, functions of design 
variables w. A typical function of displacement (e.g., a constraint) will be represented as 


g = g(U,v), U = U(y) (15) 

Analytical calculations of derivatives of displacements and their functions have been per- 
formed by three methods: the direct or design space method, the atjoint variable or state 
space method, and the virtual load method. The virtual load method is a special case of 
the direct method. Both the direct and adjoint methods begin with the differentiation of 
Eqs. (14) and (15). 


dv ov ov 


( 16 ) 
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d U 
dv 


( 17 ) 


d£ 

du 




The direct method is to solve Eq. (16) for dU /dv and substitute dU Jdv into Eq. (17). 
Equation (16) needs to be solved once for each design variable v so that the direct method 
is costly when the number of design variables is large. 

The adjoint variable or state space method starts by defining a vector of adjoint 
variables that satisfies the equation 



(18) 


where dg/dU is sometimes referred to as the dummy load vector. (If g is a particular 
displacement component, then dg/dU corresponds to a force of unit magnitude in the 
direction of the component.) Then using Eqs. (16 18), 


dp 

du 


= 2i + 

dv +X 


R v 


(19) 


The adjoint variable method requires the solution of Eq. (18) once for each function g. 
Therefore, if the number of functions is smaller than the number of design variables, the 
adjoint variable method is more efficient and, conversely, if the number of design variables 
is smaller, the direct approach is more efficient. Both the direct and adjoint methods 
involve fewer computations than the finite difference approach, which requires repeated 
factorization of the stiffness matrix, whereas the direct and adjoint methods require a 
single factorization with several right-hand sides. 


Calculation of dK/dv 

An important computational task in the adjoint and direct methods is the calculation 
of dK/dv. If the structural model contains only elements whose stiffness matrix is pro- 
portional to v (such as rods where v is the cross-sectional area or membranes and shear 
panels where v is the thickness), dK/dv is a constant matrix. But for elements having 
bending stiffness such as beams and plates, the stiffness matrix is a nonlinear function of 
the cross-sectional dimensions, and the stiffness matrix derivatives are not easily evaluated. 
The difficulties associated with shape design variables are even more severe. Analytical 
expresssions for derivatives of the stiffness matrix are cumbersome and more expensive 
to evaluate than the stiffness matrix itself. Furthermore, coding analytical derivatives of 
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stiffness matrices with respect to all possible design variables is a formidable task, espe- 
cially that in many cases users of structural analysis software that does not have sensitivity 
capabilities do not have access to the source code of the software. For these reasons, the 
preferred approach by most analysts is to compute dK /dv by finite differences. This com- 
bination of analytical derivative experssions such as Eq. (19) coupled with finite-difference 
evaluation of the stiffness matrix is known as the semi-analytical method. 

Unfortunately, the semi-analytical method is prone to large errors for some shape de- 
sign variables. The problem was explained in Ref. 20 by noting that Eq. (16) treats the 
sensitivity of the displacement vector as the solution of a structural analysis problem with 
the load replaced by the pseudo-load vector R v . This presupposes that the derivative of 
the displacement vector is a legitimate displacement vector itself, which is not always the 
case. A simple example when the derivative of the displacement vector is not a legitimate 
displacement is a nearly incompressible material (Poisson’s ratio close to 0.5). The deriva- 
tive of the displacement with respect to shape changes, treated as a displacement field, 
would typically represent large volume changes. Thus the pseudo-load vector, R v , would 
need to have extremely large components to extract such large volume changes from a 
nearly incompressible material. In such a case, the small truncation errors associated with 
the finite-difference calculation of the pseudo load are greatly amplified with a resulting 
very poor accuracy of the semi-analytical sensitivities. A similar phenomenon can occur 
for shape changes in bending problems, such as those associated with beams, plates and 
shells. The sensitivity field is often dominated by shear deformations. Since it is very 
difficult to force a slender beam or a thin shell to undergo large shearing deformations, we 
again require very large pseudo loads with disastrous effects of small errors in these loads. 

Calculation of Second Derivatives 

Second derivatives of displacement and constraint functions are used for approximate 
analysis, and for the calculation of derivatives of optimal solutions. Such derivatives may 
be obtained by differentiating Eqs. (16) and (17), for example, 


d 2 U dRv dRv d U 
K dv 2 dv + dU du 

( 20 ) 

<p 9 _ 9 ( d 2± Y *£ , (*l\ t 

dv 2 dv + \dUdv) d v \dU ) dv 2 
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However, for m design variables there are m(m -f l)/2 second derivatives, and Eqs. (20) 
need to be solved for that many right-hand sides. It is possible to proceed with a more 
efficient approach to use Eq. (18) to obtain 


= &9_ ( £9 \ df/ r f9R v dR v dU\ 

dv 2 dv 2 \dUdv) dv V dv dU dv ) 


( 21 ) 


This approach requires the solution of Eq. (16) for all the first derivatives and Eq. (6) for 
all vectors of adjoint variables. 

Stress Derivatives 

The stresses in an element may be obtained from the displacements using 

o = SU — GT (22) 

where a is a vector of element stresses, T is an element temperature, and S and G are 
stress-displacement and stress- temperature matrices, respectively. 

Derivatives of stresses may be obtained by differentiating Eq. (22). 


— = S— + —U- —T 
dv dv dv dv 


(23) 


For finite elements such as rods, membranes, and shear panels, S and G are independent 
of v, and stress derivatives are obtained by simply substituting dU/dv for U and T = 0 
in Eq. (22). For bending-type elements, S and G may be functions of v and the complete 
expression must be used; see Camarda and Adelman. 21 


Derivatives of Nonlinear Response 

In the case of nonlinear analysis, the equations of equilibrium may be written as 

P(U,v) = fiF(v) (24) 

where P is the internal force generated by the deformation of the structure, and 11 F is 
the external applied load. The load scaling factor /i is typically used in nonlinear analysis 
procedures for tracking the evolution of the solution as the load is increased. This is useful 
because the equations of equilibrium may have several solutions for the same applied loads. 
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By increasing fj. gradually we make sure that we obtain the solution that corresponds to 
the structure being loaded from zero. 

Differentiating Eq. (24) with respect to the design variable v we obtain 


dU_ _ dF_dF_ 
du ^ dv dv 


where J is the Jacobian of P at U, 


(25) 


Jki = 


m 

dUi 


(26) 


often called the tangential stiffness matrix. 

The direct method for obtaining dg/dv is to solve Eq. (25) for dU/dv and substitute 
into Eq. (17). The matrix J is often available from the solution of the equations of 
equilibrium when these are solved by using Newton’s method. Newton’s method is based 
on a linear approximation of the equations of equilibrium about a trial solution U 

P(U, v ) + J(U, v)(U - tf) « (27) 

Equation (27) solved for U, typically provides a better approximation to U than U. This 
new approximation replaces U in Eq. (27) for the next iteration, either with an updated 
value of J (Newton’s method) or with the old value ( modified Newton’s method). The 
iteration continues until convergence to a desired accuracy is achieved. If the last iterate 
U, for which J was calculated, is close enough to U, then that J can be used for calculating 
the derivative of U. 

The adjoint approach is very similar to that used in the linear case. The adjoint vector 
A is the solution of the equation 


J T A = 


dg_ 

dU 


Then it is easy to check that we obtain 


= + \ 77 lf — - —) 

dn dv K ^dv dv } 


Sensitivity of Eigenvalues and Eigenvectors 


(28) 


( 29 ) 
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Distinct Eigenvalues 


The general problem is to compute derivatives of eigenvalues and eigenvectors with 
respect to design variables or system parameters. For reference purposes, the most general 
case considered is the following eigenvalue problem: 

AX = A BX (30) 

Y t A = \Y t B (31) 

Y t BX = 1 (32) 

where A is an eigenvalue (generally complex). The generally nonsymmetric real n x n 
matrices A and B are assumed to be explicit functions of a set of design variables v, 
and X and Y axe right and left eigenvectors, respectively. The first result on eigenvalue 
derivatives was published by Jacobi, 22 who developed the result for the special case of 
symmetric A, and B = I 


— = Y t — X 
dv dv 

Wittrick 23 applied Jacobi’s formula for the case of a symmetric matrix to the derivatives of 
buckling eigenvalues and presented results for the change in buckling loads of plates with 
respect to aspect ratio and thickness. 

Fox and Kapoor 24 and Fox 25 considered the special case of symmetric A and B ma- 
trices. For eigenvalues their formula is 



« W"-a£) 

dv \dv dv ) 


X 


(34) 


in which it is assumed that the eigenvectors are normalized such that 

X t BX = 1 (35) 

For eigenvector de rivatives, two methods are presented by Fox and Kapoor. The first is 
to differentiate Eq. (30), giving a set of simultaneous equations for the eigenvalue and 
eigenvector derivatives. Differentiating the eigenvalue problem of Eq. (30) gives 




- ^ B -\^f)x 

\dv dv dv J 


( 36 ) 
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The matrix A — A B is singular since A is an eigenvalue. The set is solvable only after alge- 
braic manipulation, which destroys the banded nature of the equations, a point that arises 
later in connection with another method. The second method for eigenvector derivatives, 
developed by Fox and Kapoor, is to expand the derivative as a series of eigenvectors. Thus, 
for the zth eigenvector 


dXj 

dv 


71 

fc=l 


(37) 


The coefficients aik are obtained by substituting Eq. (37) into Eq. (36). In principle, 
it is necessary to use all n modes in the expansion of Eq. (37). However, as with the 
modal method generally, it should be possible to obtain reasonable results with fewer than 
n eigenvectors. A modification of the method of Fox and Kapoor which has exhibited 
faster convergence 26 is denoted the modified modal method 27 . This method represents 
the eigenvector derivative as 


dXj 

dv 



+ 'Yl ^ikXk 


1 


(38) 


where 


is denoted a “psuedo static” solution which satisfies the equation 



d\ dA 

dv dv 



Xi 


(39) 


The coefficients a,* are obtained by substituting eq. (38) into eq. (36). 

Rogers 28 and Stewart 29 derived sensitivity formulas for eigenvalues and eigenvectors 
of the general problem [Eqs. (30) and (31)]. For eigenvalues the equation is 


d A 
dv 


Y t 


(dA _ dB\ 
V dv dv J 


X 


(40) 


Rogers expressed the eigenvector derivatives as an expansion in terms of the eigenvectors 


dXj 

dv 


fc=l 


dY 

dv 


E 6 <‘ y ‘ 


fc=i 
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(41) 


The coefficients Qifc &nd 6^ axe computed by substituting Eqs. (41) into an expression 
obtained by differentiating the eigenvalue problem and combining it with appropriate or- 
thogonality conditions. 

An alternate method for calculation of eigenvector derivatives for the symmetric prob- 
lem is due to Nelson. The method of Nelson is to represent the eigenvector derivative 
as 

^ =v+cx («> 
where V is the solution of a reduced version of Eq. (36) obtained by deleting the fcth row 
and column from A — XB (where k is chosen to correspond to the maximum component 
of X ) and setting the fcth component of V equal to zero. The multiplier c is evaluated 
by substituting Eq. (42) into Eq. (36). This method has certain advantages over previ- 
ous eigenvector derivative techniques: it requires only the eigenvalue and eigenvector for 
the mode being differentiated, and the equation for V retains the banded character of 
coefficient matrix unlike the algebraic methods (e.g., Fox and Kapoor). 

Repeated Eigenvalues 

The sensitivity of repeated eigenvalues has been a focus of recent interest, even though 
the eigenvalues arc not differentiable and only directional derivatives can be found. For the 
real symmetric case, a generalization of Nelson’s method which preserves the bandedness of 
the matrix was obtained by Ojalvo 31 and amended by Mills-Curran 32 and Daily 33 . These 
methods compute the derivatives of the m eigenvectors corresponding to eigenvalues of 
multiplicity m. As stated by Dailey (see also Lancaster 34 ), when the eigenvalues are 
repeated and a design variable is perturbed, the eigenvectors “split” into as many as m 
distinct eigenvectors. We seek the derivatives of these distinct eigenvectors which “appear” 
with design variable perturbation. Using Dailey’s notation, define the eigenvalue problem 

KX - MX A, (43) 

where X contains the m eigenvectors cited previously, and 

A = A/, (44) 

where A is the repeated eigenvalue and / is the identity matrix of order m. The normal- 
ization condition, Eq. (35) is now 


X t MX = /. 


(45) 
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The eigenvectors which appear as a result of the splitting are contained in a matrix 
denoted Z which is related to X as follows 

Z = * (46) 

where is a set of orthogonal vectors to be determined. To simplify the notation we 
Consider a single design variable v, and denote derivatives with respect to that design 
variable by a prime. The technique for calculating Z' as contained in Daily is outlined 
next. The vector and the derivative of the multiple eigenvalues A' are obtained as 
solutions of the following eigenvalue problem 

D = A', (47) 

where 

D = X t {K' - \M')X, (48) 

with a normalization condition 

T = /. (49) 

Next in a manner analogous to Nelson 30 let 

Z' = V + ZC, (50) 

where V is the solution to 

(K - A M)V = (AM' - K')Z + MZA\ (51) 

(numerically obtained by removing m rows and columns from K - AM using the strategy 
described in Reference 32) and C is a matrix which is obtained as the solution to the 
equation 


CA'-A'C+^A" = -Z T (K'- AM') V-Z t (M'Z+ MV) A' + ^ Z T (K" - AM") Z = R. (52) 

Equation (52), which requires substantial algebraic manipulations for its derivation, deter- 
mines the matrix C and the matrix of second derivatives of the eigenvalues A". Fortunately 
A" is diagonal and CM — MC always has zero on the diagonal. Therefore, we can solve 
for the matrix C separate from A", and the latter matrix only needs to be calculated if it 
is needed for some other purpose. 
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Using Eqs. (45), (46) and (49) we have 


Z T MZ = t X t MX 

= I. 

(53) 

Differentiate Eq. (53) and use Eq. (50) to obtain 



C + C T = -V t MZ - Z t MV - 

z t m'z = g, 

(54) 

from which 

= 2 Qii 

The non-diagonal elements of C are 


(55) 

^ " v_A' 

a; ^ A'. 

(56) 


For the case where A' = A' i ^ j, eq. (56) may not be used. The situation here is 
that the eigenvalues are not “splitting” when the design variable is perturbed because the 
design variable is affecting both in exactly the same way. In such a case, Z' is not unique 
and any values of and Cji satisfying Eq. (54) may be used. Dailey proposes the choice 
Cij = Cji = whenever A' = A'-. 

Before leaving the topic of derivatives associated with repeated eigenvalues, we note 
the limited utility of such derivatives. For example, the eigenproblem is differentiable in 
terms of a single parameter, but not as a function of several. This may be demonstrated 
by the example where the matrix 


A = 


2 + y 

x 


X 

2 ' 


(57) 


The eigenvalues of A are 

Ai ,2 = 2 + y/2 ± v^x 2 + y 2 /4. (58) 

At x = y = 0, the eigenvalues are repeated and dX/dx = ±1, dX/dy = 0, 1. However, the 
eigenvalues arc not differentiable as a function of both x and y, that is the relation 


dX dX 

AX = ai Ax + ai Ay 


( 59 ) 
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does not hold. Therefore, the utility of the partial derivatives is questionable. The eigen- 
vectors are also discontinuous at (0,0). This can be checked by noting that at (e, 0), the 
eigenvectors are (1,0) and (0,1) and at (0,e) they are (1,1) and (1,-1) no mater how 
small e is. 

Sensitivity derivatives for nonlinear eigenvalue problems 

In flutter and nonlinear vibration problems, we encounter eigenvalue problems where 
the dependence on the eigenvalue is not linear. For example, Bindolino and Mantegazza 35 
consider aeroelastic response problem which produces a transcendental eigenvalue problem 
of the form 

A(A,v)AT = 0. (60) 

Differentiating Eq. (60) we get 


dX 


dXdA 
+ dvdX 


- = -**X. 


(61) 


dv dv d\ dv 

Using the normalizing condition X m = 1 we can solve Eq. (61) for dX/dv and dX/dv. 
Instead, it is also possible to use the adjoint method, employing the left eigenvector Y 
satisfying 

Y T A = 0, Y m = 1 (62) 

we obtain 


dX 

dv 


Y T ^X 

dv 

itX ■ 

y T ™X 

dX 


(63) 


A common treatment of flutter problems is to have two real parameters representing 
the approach of the frequency and speed as an eigenpair instead of one complex eigenvalue. 
For example in Murthy 36 , Eq. (60) is replaced by 

A(M t u>,v)X = 0, (64) 

where the Mach number, M, and the frequency, u>, are real parameters. Using this ap- 
proach, differentiate Eq. (64) and premultiply by Y T to get 


/m 1ST + / ‘ 


du) 

“ Th = ~ fv ' 


(65) 
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where 


dA 


dA 


.dA 


fM = y m x ’ ^ = y /• = 


c?u; 


dv 


( 66 ) 


Multiplying Eq. (65) by f u (the complex conjugate of /„,) we get 


, f dM . , l2 duj T 

Im/u + \ fu \ q v — fwfv 


(67) 


The second term in Eq. (67) as well as are real, so by taking the imaginary part of 
Eq. (67) we get 


dM 

dv 


Im(Ufv) 

'(Y T %X) (yT!£x ) 

im 

\(Y T §a x ) {y t $x) 


Next, multiplying Eq. (65) by /m and following a similar procedure gives 


du) 

dv 


Im 


(Y^X) (Y T j$x) 


Im 


(Y T m x ) (? T S£*) 


( 68 ) 


(69) 


Rudisill and Bhatia 37 have a derivation of the flutter eigenpair that employs the 
reduced frequency and flutter speed as the eigenpair and provides also second derivatives. 

It is possible to treat in a similar manner the case where the nonlinearity is in X 
instead of in A. For example, Hou et al. 38 treated the nonlinear vibration problem 

K(X)X - A MX = 0 (70) 

Differentiating Eq. (70) with respect to v we obtain 


fT \m\®X d\.. Y f dK dM 

( J — -q^MX = — X-g^]X 

where J is the tangent stiffness matrix whose components are given as 


(71) 


jij — Kij + 

k 


dKik 

dXj 


X k 


(72) 
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Equation (71) can now be solved for eigenvector deriviatives using Nelson’s method. For 
eigenvalue of derivatives use the left eigenvector satisfying 

Y t (J-XM) = 0 

Premultiply Eq. (71) by Y T to obtain 

d\ Y T (%-\e$)X 

dv Y T MX 

Sensitivity of Transient Response 

General 

The discussion of sensitivity analysis of transient structural response is usually based 
on the equations of motion written as a system of second-order differential equations. 
However, this form obscures the similarity of structural sensitivity analysis to sensitivity 
analysis in other fields where first -order differential equations are employed and is also less 
compact than a first-order formulation. For these reasons the initial discussion will focus 
on a system of first-order differential equations, and the case of a second-order system will 
be limited to linear structural dynamics, following the general discussion. We start with a 
system of first-order ordinary differential equations of the form 


(73) 


(74) 


U = F(U,t,v) 

U(0) = U 0 


(74) 


where U is the response, F a vector of functions, t time, and v a typical design parameter; 
a dot denotes differentiation with respect to time. In many structural applications the 
left-hand side of Eqs. (56) is A U, where A is a matrix, and the methods discussed below 
are also applicable to that more general form. 

Direct Method 

The direct method of obtaining sensitivity derivatives is based on differentiating 
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Eqs. (74) to obtain 


d U_ _ jdU_ _ dF_ 
dv dv dv 

( 75 ) 


d U 
dv 


( 0 ) = 0 


where the Jacobian J is df jdU . Note that Eqs. (75) are a system of linear differential 
equations, even if the original system, Eqs. (79), is nonlinear. Often derivatives of the 
entire vector U are not required. Instead it is necessary to obtain the derivatives of a 
function of U of the form. 

f tf 

g(U,v)= p(U,t,v)dt (76) 

Jo 


where p is a functional representation of a time-dependent constraint and tj is a final time 
for the response calculation. The direct approach obtains dg/dv as 


d g = [*' \dp + (dp_\ T dU 

dv J 0 dv \dU J dv 


(77) 


where dU/dv is calculated in Eqs. (75). 

Green’s Function Method 

Equation (75) have to be solved once for each design variable and are costly when 
the number of design variables is large. When the number of design variables is larger 
than the dimensionality of U, then the Green’s function approach 9 is more efficient than 
the direct approach. An application of this approach is sensitivity analysis of transient 
structural response, when the response is computed using reduction techniques such as 
modal analysis. The sensitivity derivatives dU/dv is written as 


dU 

dv 


(0 



(78) 
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where the Green’s function K satisfies (recall that the dot denotes d/d£) 


K(t,r) = 0, t<r 

K(t,t) = I (79) 

K(t,r) — J(t)K(t,r) = 0, t>T 


The efficiency of the Green’s function approach is partly governed by the method used to 
integrate Eqs. (79). A large amount of work on the efficient implementation of the Green’s 
function approach has been performed by Rabitz and co-workers. 39 


Adjoint Variable Method 

Further improvements in efficiency may be possible if less information is needed. If 
instead of the derivatives of the entire vector U, only those of a few functionals [e.g., 
Eq. (76)] are required, then an adjoint variable method is called for. The adjoint variable 
approach solves first for the adjoint vector A from the differential equation 


A + J t A 


dp 

~dU 


A (if) = 0 


(80) 


It is shown by Haftka, Giirdal and Kamat 40 that 


dg 

dv 




(81) 


Equations (80) are a set of linear differential equations that is integrated backward from 
tf to zero. As in the steady-state case, the adjoint variable approach is preferred over th e 
direct approach when the number of functionals is less than the number of design variables. 

Finite Difference Method 

For sensitivity analysis of static response, the finite difference approach is almost al- 
ways inferior to analytical methods. For the calculation of derivatives of transient response, 
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this is not always the case. When explicit methods are used for integrating the differential 
equations, the linearity of the sensitivity equations does not constitute a computational 
advantage. Therefore, for the case of explicit integration, the finite difference approach is 
often computationally superior to the direct method (see Ref. 41). When implicit integra- 
tion t echni ques are used, the finite difference approach is less attractive computationally 
but remains easier to implement than the direct approach. 

Linear Structural Dynamics 

For the case of linear structural dynamics it may be advantageous to retain the second- 
order equations of motion rather than reduce them to a set of first-order equations. It is 
alsn common to use modal reduction for this case. In this section we discuss the application 
of the direct and adjoint methods to this special case. The equations of motion are written 
as 

MU + CU + KU = F(t) (82) 

Most often the problem is reduced in size by expressing U in terms of m basis functions 
^ i — x, . . . m where m is usually much less than the number of degrees of freedom of the 
original system Eq.(82) 

U = $Q (83) 

where $ is a matrix with as columns. Then a reduced set of equations can be written 


as 

MQ + CQ + KQ = ?(t) 

(84) 

where 


M = = K = * t K$, F = $ t F 

(85) 


When the basis functions are the first m natural vibration modes of the structure scaled 
to unit modal masses $ satisfies the equation 


= 0 ( 86 ) 

where a is a diagonal matrix with the ith natural frequency c* in the ith row. In that case 
K = Q 2 and M = I are diagonal matrices. For special forms of damping, the damping 
matrix C is also diagonal so that the system Eq. (84) is uncoupled. After Q is calculated 
from Eq. (84) we can use Eq. (83) to calculate F. This method is known as the mode- 

displacement method. : 7 . . -- 

When the load F has spatial discontinuities the convergence of the modal approxi- 
mation, Eq. (83), can be very slow. The convergence can be dramatically accelerated by 
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using the mode acceleration method. The mode acceleration method can be derived by 
rewriting Eq. (82) as 

U = K~ l F - K~ l CU - K~ l MU (87) 

The first term in Eq. (87) is called the quasi-static solution because it represents the 
response of the structure if the loads are applied very slowly. The second term and third 
terms are approximated in terms of the modal solution. It can be shown (e.g., Greene 42 ) 
that K~ l can be approximated as 


j 




K~ l = ( 88 ) 

Using this approximation for the second and third terms of Eq. (87) we get 

U fvK-'F -m- 2 CQ-m~ 2 Q (89) 

This approximation is exact when $ contains the full set of vibration modes. Note that Q 
and Q in Eq. (89) are obtained from the mode-displacement solution, Eq. (84). Therefore, 
there is no difference in velocities and accelerations between the mode-displacement and 
the mode acceleration. 

In considering the calculation of sensitivities we treat first the mode-displacement 
method. The direct method of calculating the response sensitivity is obtained by differen- 
tiating Eq. (84) to obtain 


where 


m^ + c^ + k^- = r 

dv dv dv 


„ dF dM - dM - dK 

R= Q Q Q 

dv dv dv dv 


(90) 

(91) 


The first step in forming this equation is the calculation of the derivatives of F, M, C, 
and K with respect to v. Differentiation of K yields 


dK T dK d$ T __ T dfy 

—— = 4> - 7 -$ H — ; — K$ + 4> K - — 
dv dv dv dv 


(92) 


with similar expressions for the derivatives of M, C, and F. The calculation is simplified 
considerably by using a fixed set of basis functions $ or neglecting the effect of the change 
in the modes. In many cases the error associated with neglecting the effect of changing 
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modes is small. When this error is unacceptable we have to face the costly calculation of the 
derivatives of the modes needed for calculating the derivatives of the reduced matrices, such 
as Eq. Fortunately it was found by Greene 42 that the cost of calculating the derivatives of 
the modes can be substantially reduced by using the modified modal method 26 , keeping 
only the first term in this equation. This approximation to the derivatives of the modes 
may not always be accurate, but it appears to be sufficient for calculating the sensitivity 
of the dynamic response. 

For the adjoint method we consider a function in the form of Eq. (76) 


a{Q,v) 



p(Q,v , t)dt 


(93) 


so that 


dg 

dv 



dp dQ 
dQ dv 


)dt 


(94) 


To avoid the calculation of dQ/dv we use an adjoint vector A, and start by multiplying 
Eq. (90) by X T and integrating 




A T Rdt 


Integrating by parts we get 


(95) 


\ t m ^ 

dv 


-A 

dv 




+ A 


TgdQ 

dv 


+ 


J \x T M-X T C+X T K)^ = J ’ X T Rdt (96) 


Assuming that the initial conditions do not depend on the design variable v, Eq. (96) 
suggests the following definition for A 


MA-CA+ KX = (||) 7 \ 


A (*/) = A (t f ) = 0 


(97) 
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and then Eq. (94) becomes 



( 98 ) 


For the mode-acceleration method we consider only the direct method. We start by 
differentiating Eq. (82) and rearranging it as 


dv 


dF 

dv 


dv dv dv dv dv 


(99) 


Next wc use Eq. (89) to approximate the second term, and the modal expansion Eq. (83) 
to approximate the other terms to get 


dU 

di> 


K 


-l 


dF dK 

*r - * ri (*' F - ™->cq - 


dQ dC - dQ dM - 

Gv- — — - M<b~ $0 

dv dv dv dv v 


( 100 ) 


Finally we use the modal approximation to A' -1 , Eq. (88) to obtain 


d U 
di> 




-l 




<MT 2 $ t 


~‘hn- 2 CQ - ^#0 - Ci>~- 

dv dv dv 


+ 


( 101 ) 


K~ l 



dM 

dv 


$ 


Q-$n 


-2 


dQ 

dv 


Note that the calculation involves the solution of Eqs. (84) and (90) for Q and dQ /dv, 
followed by Eq. (101) for retrieving the dU/dv. Additional details can be found in 42 . 
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